Metastable states and T = hysteresis in the random-field Ising model on random 

graphs. 
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in I 

, We study the ferromagnetic random-field Ising model on random graphs of fixed connectivity 

(Bethe lattice) in the presence of an external magnetic field H. We compute the number of single- 
' spin-flip stable configurations with a given magnetization m and study the connection between the 

distribution of these metastable states in the H — m plane (focusing on the region where the number 
is exponentially large) and the shape of the saturation hysteresis loop obtained by cycling the field 
between — oo and +oo at T = 0. The annealed complexity £a(wi, H) is calculated for z = 2,3,4 
and the quenched complexity T.Q(m,H) for z = 2. We prove explicitly for z = 2 that the contour 
lO ■ ^Q(m,H) = coincides with the saturation loop. On the other hand, we show that £a(ti, H) is 

irrelevant for describing, even qualitatively, the observable hysteresis properties of the system. 



PACS numbers: 75.10.Nr, 75.60.Ej, 05.50.+q 

I. INTRODUCTION 
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The dynamical properties of random and frustrated systems are dominated at low temperature by the presence 
of a multitude of metastable states. This complex energy landscape is responsible for the strong hysteresis effects 
observed in a variety of physical systems (magnets, superconductors, fluids in porous media,. . . ) under the action of an 
' external driving field. Much insight into such behavior has been gained recently through studying the zero-temperature 
random-field Ising model (RFIM) in an external magnetic field, which is a simple prototype of a disordered system 
far from equilibrium. The energy landscape of the model may be explored by changing the field adiabatically, which 
O \ results in a series of discrete Barkhausen jumps (avalanches) in the magnetization, each jump corresponding to a 
sudden move to a new metastable configuration. A remarkable feature of the RFIM at T = is the existence of a 
| disorder-induced phase transition in the magnetization hysteresis loop obtained by cycling the field from one saturated 
• state to the oppositely magnetized state and back. In three and higher dimensions, there is a transition between a 
continuous loop at strong disorder and a discontinuous loop with a macroscopic jump in the magnetization at low 
disorder^. The jump (corresponding to a macroscopic avalanche in the system) disappears at some critical value of 
the disorder where universal scaling behavior is observed. This noncquilibrium phase transition has been studied in 
detail by renormalization group and by numerical simulations on hypercubic lattices . Exact analytical calculations 
have also been performed in one dimension^ and generalized to a Bethe lattice^]: in this case, the transition only 
exists for z > 4, where z is the connectivity of the lattice. Such a transition has been observed in thin Co/CoO 
films|f| and Cu-Al-Mn alloysQ, and we have recently suggested^ that it may also be associated to the change in the 
adsorption behavior of 4 He in dilute silica aerogels 8J. In all these systems, activated processses are too slow to allow 
thermal equilibration on experimental time scales, which makes the T = RFIM a relevant model for understanding 
I ' the out-of-equilibrium behavior. 

There is an expected connection between the shape of the T = hysteresis loop and the distribution of the 
metastable states in the field-magnetization planejij. This issue, however, has only attracted little attention so 
farpjl EI. For the RFIM with ferromagnetic interactions, the probability to find metastable states outside the 
hysteresis loop vanishes in the thermodynamic limit, as shown in Appendix A. But what is the situation inside the 
loop ? Is the entropy density (associated with the number of metastable states) positive everywhere ? The situation 
is especially intriguing in the low-disorder regime where a region around the origin becomes inaccessible to any field 
history that starts from the saturated states [lj- The aim of this work is thus to calculate the number of metastable 
configurations with a fixed magnetization m per spin as a function of the external field H . We consider a Gaussian 
distribution of the random fields whose width R measures the "strength" of the disorder, and, in order to analyze 
the problem analytically (at least partially), we work on a Bethe lattice defined here as a random graph of fixed 
connectivity z |13j, i.e. a random collection of ./V vertices, each vertex being connected to z other vertices. Such a 
graph has the structure of a branching tree locally and small loops are rare, with typical size of order log N. Spin 
models on random graphs (both with fixed and fluctuating connectivity) have attracted much attention in recent 
years because they can be solved exactly though each spin interacts with a finite number of othe r sp ins like in real 
short-range, finite-dimensional systems. The metastable configurations of the Ising ferromagnet [141 Il5l Ha | and of spin 
glass models^! El have been studied in this framework, but, to the best of our knowledge, the ferromagnetic RFIM 
has not yet been considered (except for 2 = 2 and in zero external field p^K Since there is no plaquette frustration 
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in this case (contrary to spin glasses), loops should not play any significant role in the large- N limit and we expect 
the number of metastable states on a t ypi cal random graph with connectivity z to be the same as in the interior of a 
Cayley tree with branching ratio z — 1 • 

As is usual with disordered systems, one has to distinguish between annealed and quenched averages, the disorder 
average being here over both the random-graph and the random-field distributions. The quenched average is the 
appropriate quantity describing the behavior of the typical disorder realizations. Specifically, we want to compute 
both the average number Af(m, H) and the typical number e ln ^( m ^ H ) f single-spin-flip stable states with a given 
magnetization m. Since these numbers are expected to scale exponentially with the number N of spins (at least in 
some region of the H — m plane), the basic quantities are the associated entropy densities, i.e. the annealed complexity 
Y,A(m,H) — limAT^oo 1/A ln7V(rn, H) and the quenched complexity T,q(iti,H) = limjv^oo 1/A h\N{m,H). In the 
following, we present the results for the annealed complexity on random graphs with connectivity z = 2,3 and 4. 
We shall see that the qualitative behavior changes for z > 3, although the recursion equations for the RFIM on 
a Bethe lattice predict that the macroscopic jump in the hysteresis loop only occurs for z > 40. This illustrates 
that the annealed quantity is unable to describe, even qualitatively, the typical behavior of the system: indeed, the 
probabilistic calculation of Ref . Q describes the behavior of the magnetization which is a self-averaging quantity that 
converges to its typical value when N — + oo. (Note that the crucial distinction between annealed and quenched 
averages was neglected in Ref.|ll]|.) The analytical computation of the quenched complexity requires the introduction 
of replicas and is thus more complicated (alternatively, one could also use the cavity method as in Ref. .16:]). For 
this reason, we carry out the complete calculation of T,Q(m,H) for z = 2 only. We then show explicitly that the 
contour £q(to, if) = identifies with the saturation loop of the one-dimensional RFIM previously calculated by a 
probabilistic method 3] : the typical number of metastable states thus ceases to be exponential with the system size 
exactly on the saturation hysteresis loop. To our knowledge, this is the first time that such a result is proven exactly. 

The rest of the paper is organized as follows. In Sec. II we introduce the model and describe the averaging 
procedure over the random graphs. In Sec. Ill we present the general formalism for computing the annealed and 
quenched complexities on such graphs. The analytical and numerical results for the one-dimensional chain (z — 2) are 
discussed in Sec. IV. In Sec. V we present the results for the annealed complexity for z — 3 and z — 4. We conclude 
in Sec. VI by discussing the possible behavior of the quenched complexity for z > 4 in the low-disorder regime. The 
absence of metastable states outside the saturation hysteresis loop is proven in Appendix A. Appendices B and C are 
devoted to technical details of the quenched calculation for z = 2. 

II. MODEL 

We study the ferromagnetic RFIM described by the Hamiltonian 

H = -J^riijSiSj -Yj( H + h)si , (1) 

i<j i 

where J > 0, Si (i = 1, . . . , N) are Ising spins placed on the vertices of a random graph, and the random fields hi are 
drawn from the Gaussian probability distribution P(hi) — exp(— hf /2R 2 ) / ^/2ttR. The variable is equal to 1 if the 
sites i,j are connected and to otherwise. 

Following Ref.|l7|. the random graphs are built as follows. First, the elements mj of the connectivity matrix (with 
fiij — riji, i < j) are chosen at random from the probability distribution Viriij) = (1 — 7 / 'N)6(n%j', 0) + / N)8{ni 3 ■; 1) 
where 7 is some arbitrary number of order one and 5{p\ q) is the Kronecker-<5: this defines an ensemble of graphs where 
the local connectivity is a Poisson variable in the large- N limit, like in the Viana-Bray spin glass model[2(|. One 
then selects the subset of graphs where the connectivity of each site is equal to z. The average value of a dynamical 
variable A is thus given by 

1 N 

A= M(iV,z,7) ( ^n^E^^))7)^ (2) 

where (.) 7 and (.)/, denote the averages with respect to the distributions V{riij) and P(hi), respectively. M{N, 2,7) 
is the average number of graphs with connectivity z that are generated for a given 7. This quantity was calculated in 
Ref.fnj in the large- A'' limit: 

lim — lnM(N,z,j) = -(lnz + ln 7 - 1) -lnz! - - . (3) 
N— >oo A 2 2 



3 



By construction, the value of A is independent of the choice of 7 in the thermodynamic limit N — > 00. In the following, 
we drop the explicit dependence on z and 7. 

III. FORMALISM 

In this section we present the general formalism for computing the number of metastable states of the Gaussian 
RFIM on random graphs. Since the quenched calculation is ac hieved by in troducing replicas and using the identity 
lux = lmi n ^o(x n — l)/n , we compute directly the average Af(m,H) n for a general integer n, and then treat 
separately the annealed (n = 1) and quenched (n — » 0) cases. Our method is inspired from Ref.[l5[ where a global 
order-parameter function is introduced that allows to adress the quenched calculation. 

A. Basic equations 

A configuration is said to be metastable when its energy cannot be decreased by flipping a single spin, a definition 
consistent with the T = Glauber dynamics used in the previous studies of the nonequilibrium RFIM [J, La-Ll Such 
a configuration is a local minimum of the energy surface where each spin i is aligned with its local field fi defined by 

fi = J n ii s i + h i + H ■ ( 4 ) 
j 

The nth moment of the number of metastable states averaged over disorder is thus given by 



• V < H >" - -iW (5) 



with 

N 



Si— ±1 i—1 j^i 

= «Tr H[e(f?sf)]6(J2nif,z)h)h (6) 

i,a j^i 

where Tr stands for the summation over all possible values of the n-replicated spins sf (a = 1, n), /? = J J^j n ij s j + 
hi + H, and Q(x) is the Heaviside step function (i.e. Q(x) — 1 if x > and if x < 0). Since we consider a continuous 
distribution of the random fields, there is no need to specify the value 0(0): the local field /j is almost surely different 
from zero and the system has no marginally stable states whose energy does not change under any single-spin flip. 
We then use integral representations of the Dirac-5 and Kronecker-J to rewrite Eq.(6) as 

V n (N,H) = ((Tr J [dxWx* - /? )Q(x^t)] 6(£ n ^ 2 )>7>» 

i,a j^i 

= ((Tr J ft [<dyfe(<<)]e^*^-^- f ^5(X;n ij -;^)) 7 )h 



= « Tr / II [dx?d»fdA,e(x? fl ?)]eE*^-(*«-ft)-^C^« W -»)]) 7 ) h (7) 

where the bold symbols denote vectors in replica space, e.g. Xj = (x},xf, x"). The integration ranges are [—00, +00] 
for xf, [— ioo, +ioo] for t/?, [0, 2iir] for Aj, and for brevity the normalisation factors l/(2iw) are absorbed in the 
infinitesimal elements dyf and d\i. Inserting the definition of ff (with h.^ = hil,Ji = HI) and performing the 
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average over the distribution Viriij) yields 

V n (N,H) = (Tr Jll \dx a l dy a l d\ l Q{x a l s a l )y^^ z+y '-^'-^- ll)] 



nfl _ _L i _L e -Ai-A j -J(y 1 .s i +y J .s 4 )N\ 
N N " 



(Tr j H [dx?dy:d\Mx-s?) 



^JA^+y^x.-hi-H)] 



exp(_2 ? + ^E^ Al ~ A ^ J(y - s,+yj ' s,) ))^ ( g ) 



where terms that are negligible in the thermodynamic limit have been dropped in the last exponential. 

Like in Ref.0] we now introduce an "order-parameter" function c(<t,t) = (1/N) ^{sf, tr) exp[— (A, + Jy.^r)], 
where tr and r are binary vectors in replica space (i.e. a a ,r a = ±1). Then, 

ex P [-f + -Lj2^- J ^ +y ^] = «p[-^ + f £c( CTi t)c(t, CT )] , (9) 

a relation that we use to decouple the sum over the site labels i and j and perform the trace over a single spin only 
in each replica. As outlined in the Appendix of Ref.|15|. the In- variable function c(er, t) may be inserted in Eq. 
(8) via integrals over (5-functions represented by Fourier integrals over other auxiliary functions. After some formal 
manipulations, we obtain 

V n (N,H)= f [dc( CT ,T)]exp{-^£ c(CTjX)c(r , CT) _2^ } 
"* cr,T tr,T 

\J2 [ P(h)dh j dxdydAexp{Az + y.(x-h-H)+ 7 e- A ^c(x,cr)e- Jy ' r }]Je(a; a (7 a ) ** (10) 



X 

tr 



which may be computed in the thermodynamic limit as a saddle-point integral. This leads to 



lim ^V n {N,H)= ^{-J^c^TWr.ffJ-J+AH}, (11) 

N^oo W C(<T,T) 2 ~ 



tr,T 

with 



A[c]=ln^ / P(h)dh ( dxrfydAexp{Az + y.(x-h-H) + 7 e- A ^c(T,cr)e- Jy ' r }]Je(a; a cr a ) . (12) 

tr T a 

Using the identity fLU] 

J- / dAe^+^ A = £L , (13) 

we get 

A[c] = ln^ +ln^ / P(h)dh f dxdy e y ( x - H - h ) 
cr 

x ^ c(T 1 ,a)...c(T z ,cr)e- J y-^ +T2+ - T ^l[Q(x a a a ) , (14) 

Ti...Tj a 

and the integrations over y (giving (5 functions again) and then over x may be performed to finally obtain 

A[c]=ln^-+ln ^ c(T 1 ,tr)...c(T z ,tr) j P(h)dh\\Q[{H + h + J^rf)o a ] . (15) 

tr,T 1 ...T z a 1=1 

Note that the random field distribution P(h) has not yet been specified at this stage. 
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It is straightforward to generalize the preceding calculation for a fixed magnetization m. T> n (N, m, H) is now defined 



as 

N 



V n (N, m, H) = « [ II ©(^WE Sl - Nm)] "<5(^ n l3 - z)) 7 ) h (16) 

5i=±l i— 1 i 

and the <5-function constraint on . Sj may be relaxed by introducing a Lagrange multiplier g coupled to the magne- 
tization (see, e.g, Ref . |T3[ ) . The whole calculation is quite similar to the one presented above and Eq.(ll) becomes 

lim ±-lnV n {N,m,H) = max {- J V c(<x, t)c(t, tr) - J - m V 3 a + A[c, g]} , (17) 
*r-«JV c(<r,r),g 2 ^ 2 a 

with 

A[c,g] =ln ll+ln c(T 1 , t r)...c(T z ,tr) ( P{h)dh\[e aa ° a Q[{H + h + J^t?)o*] (18) 

<7,Ti...*7% ^ a i=l 

where g = (.g 1 ,^ 2 , ...,g n ). As can be seen immediately, Eqs. (11) and (15) are recovered when g = and thus Eqs. 
(17-18) may be used as the starting point of further analysis. 
Extremization with respect to c(er, r) readily gives 

,-A 7 



with 



c(tr,T) = e A — — - c(xi,cr)...c(x 2 _i,cr)i ;l z (cr,r,Xi...T 2 _i) (19) 

^ Ti...T,_i 

z— 1 

F,(<7,t,t 1 ...t,_ 1 ) = y P(h)dhl[eV a ° a e{(H + h + J{r a + Y J T ?))° a \ , (20) 



which in turn yields the normalization condition 



^c(a,T)c(T,CT) = l. (21) 



^ v a c{T U (r)...c{T Zl( T) f P(h)dhY[e 9aaa e[(H + h + J^r, a K] 



Likewise, extremization with respect to g a yields 
m = 9A[c, g]/5g° 
-AT*. 

(T,T-l...T z 

= ° a <°, t)c(t, a)/ T ) C ( T > °) ( 22 ) 

<T,T <T,T 

with a = 1, ...,n. 

Inserting Eq. (21) into Eq. (17) and using Eqs. (3) and (5), we finally obtain in the large- N limit 

7V(to, H) n ~ exp [N max {A[c,gr] - m^5 a - -(In 2; + In 7) + lnz!}] . (23) 

B. Annealed complexity 



We first consider the case n = 1 corresponding to the annealed average. <r and r become standard Ising variables 
a, t and the preceding equations simplify considerably. In particular, the function c(cr, r) can only take four distinct 
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values which we denote c ++ , c + _, c_ + , and c__. Performing the integration over the random field distribution, Eq. 
(18) becomes 

A[c, g] = In J + In £ ^, a)...c(r„ <r)~ [l + erf(^±|^*)] 

o\Ti ...r 2 — ±1 

= In ^ + In [e»C_ + / + (ft±) + «"»<_/_ (^-)] (24) 



c 



where erf(x) = (2/^/tt) J q x e * dt is the error function and the functions f±(x) are defined by 



2 n W ^ i?\/2 

n— v 7 v 

Like in Ref.flij we now change to the new variables u = c ++ /c_ + , v — c__/c + _, t = c + _/c_ + , solve the saddle 
point equation for c_ + , and substitute the result in Eq. (23). We then find that the dependence on 7 disappears 
(as it must be), and the average number of metastable states with magnetization m is given in the large- N limit by 
Af(m, H) ~ exp[E A (m, H) N] with 

T, A (m,H)= max {— ln(w 2 + iV + 2t) + ln[/+ (u)e 9 + t z f_ (w)e" 9 ] -, gm} . (26) 

u,v,t,g 2 

Solving the stationarity condition with respect to g gives 
and substituting into Eq. (26) we finally obtain 

5* 1 I jji 1 

E A (m, if) = max{ - - ln(u 2 + t 2 v 2 + 2t) + —— In fju) + ^— In fJv) 
u,v,t 2 2 2 

1 — m , 1 + m , 1 + to 1 — m . 1 — 

+ z lnt In In } . (28) 

2 2 2 2 2 ; v ' 

The saddle point equations dYiA/dv — and dY>A/dt = may be solved analytically to give t = t(v) — 
fL(v)/[zvf_(v) — v 2 f'_(v)] and u = [(1 + m)t 2 v 2 + 2mt] 1 ' 2 /[l — m} 1 / 2 (the negative branch of u leads to solutions that 
have no physical meaning). On the other hand, the remaining equation, dYiA/du — 0, must be solved numerically. 

The total complexity £™ aa: (iJ) also corresponds to the maximum of E^m, H). It is given by Eq. (26) with g = 0; 
one can associate to this maximum a value of the magnetization m^(i?), which is then obtained by solving Eq. (27). 
One still has t = t{v) but the saddle point equations for u and v must be solved numerically. For H — however, 
these equations admit the solution u — v and one has to solve the unique equation t(v) — 1. As in the case of the 
pure ferromaenet |l4j . this solution corresponds to the fixed point of the transformation u— >u, v ^ u, t —* 1/t which 
leaves £™ aa: (ff = 0) invariant (when R ^ 0, this is the only solution of the saddle point equations). 



C. Quenched complexity 

We now turn to the calculation of the quenched complexity given by 



where, as usual, the order of the limits TV — > 00 and n — > is inverted and the expression of A r (m, H) n is extended 
to non-integer n to allow the limit n — > to be taken. Some ansatz concerning the form of the solution of the saddle 
point equation (19) as n — > is now required. As alluded to in the Introduction, a major difference between the 
ferromagnetic RFIM and spin-glass models comes from the role of frustration. For the Bcthe lattice spin- glass [l3j. 
frustration (due to the presence of large loops) induces a replica symmetry breaking instability similar to the one found 
in the fully-connected Sherrington-Kirkpatrick model and corresponding to the appearance of an exponential number 
of pure states. On the other hand, for the RFIM, this phenomenon does not occur in the infinite-range model 22], nor 
on finite-connectivity lattices at the mean-field level (it may be checked for instance that the number of pure states 



Af(m,H) n - 1 



(29) 
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on the Bethe lattice is indeed non-exponential |23j ) . We thus assume that the solution of Eq. (19) remains invariant 
with respect to the permutation of the replica indices when n — ► 0. Setting g a = g in Eq. (23), we then obtain 

1 z 

Eo(m, H) = max{-mj + lim —\A*(g) Qxiz + In 7) + lnz!]} 

g n->0 n 2 

= A* 1 (g*)-mg* (30) 

where A*(g) = Aq + nA*(g) + 0(n 2 ) is the value of A at the saddle-point solution c*(er,-r) and g* is the position 
of the maximum of A* (g). The second line of this equation also uses the relation Aq = |(lnz + In 7) — lnz!, which 
is a necessary condition for the n — > limit to exist and be independent of 7. Note that Eq. (30) corresponds 
to a Legendre transformation between Eq and AJ, so that the stationnarity condition dA\/dg* = m implies that 
dT,Q(m 1 H) I dm = —g*. In fact, as will be illustrated below, it is convenient to consider that both Eq and m are 
parametrized by g (to simplify the notations, the superscript * is dropped in the following), with g going from — 00 
to +00. One then has the overall symmetry Eq(— g,— H) = Eq(<?, H), m(—g,—H) = m(g,H), and the maximum 
number of metastable states is obtained by setting g = 0. 

The main consequence of replica-symmetry is that c(<t,t) is a function of the three variables s — ^2 a= i& a , 
u = X)"=i < 7 °" rQ '1 t = ELi T ° only. This suggests, following Ref . [Ti| . to introduce an integral representation of 
c(s,u,t) so as to perform an analytic continuation of this function on the imaginary axis in the small n limit where 
s, u, and t may be treated as continuous variables. It turns out, however, that the situation is not that simple here 
because of the constraints on the values of a a and r a . To illustrate the problem, let us consider the case z — 2. Eq. 
(19) then reads 

c(a, t) = je~ A Y, <r)Fn(<r, r, n) (31) 

Ti 

with 

F 2 {(T,T,T 1 ) = e 9S I P(h)dh]]e[{H + h + J(r + T?))a a ] . (32) 

Performing the integration over the random field distribution in the four regions h < —2 J — H, —2 J — H < h < —H, 
-H < h < 2 J - H , and 2 J - H < h, we obtain 



F 2 (a, t, Tl ) = L^e-^ J] S(a a ; -1) + L-lef* J] S(a a ; 1) 

a a 

+ ^T e9S II [^( <7 °; -^(^ + S(T?;r a )5(T a ;a a )] 



o9S 



J] [5{cr a ; l)6(r?l ~r a ) + Stf; r a )6(r a -a a )] (33) 



2 

a 

where a = erf((2J + H)/Ry/2), b = ed{H/R^2), and c = erf((2J - H)/Ry/2). With this, Eq.(31) reads 
- C (<r, r) = l^e-^ J] S(a a ; -1) £ c(n, -1) + ^e«" [] ^ ^ E c ( x i - X ) 

a T 1 a T 1 

+ ^ e9S E CT ) II - 1 )^ r " - T<1 ) + ^W;r a )5(r°; O] 

T 1 a 

+ E < T ^ a ) II WW 5 ~ rQ ) + t(T?;r a )6(T a ; a a )] (34) 

Ti a 

where ±1 = (±1, ...±1). The difficulty comes from the presence of the Kronecker-5's that precludes any straightforward 
use of a Fourier representation of c(t, <t) in the limit n — > 0. Because of this, we work out the self-consistent equation 
(19) before taking the small n limit, a nontrivial (if not unsolvable) task for a generic value of the connectivity. In 
fact, only the case z = 2 seems to be simple enough to be handled analytically. Inspection of Eq. (34) suggests that 
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the solution has the form 

«/|c(«r, t)=A]J S(a a ; -1) + Be 9 * ]J S(a a ;T a ) + C J] 8(a a ; 1) 

* a a a 

+ D(<t , t) [] [6(<r a l -mr a ; 1) + S(r a ; a a )j 

a 

+ E(a, t) Yl [S(a a ; l)5( T a ; -1) + S(r a ;a a )] , (35) 

a 

which in turn, together with Eq. (21), yields a set of six coupled equations for the unknown quantities 
A, A, B,C, D((t,t), E(ct,t). It is only at this stage that the n — > limit may be taken, noticing that D(a,r) 
and E(cr, t) are functions of two variables only (say, s and t) instead of three because of the constraints on cr and r 
that are imposed by the Kronecker-<5 (it is easy to check that the third variable u = n — \s — t\) . We then introduce 
the integral representations 

D(s,t) = / dxdyD(x,y)e xs+yt , E(s,t) = [ dxdyE{x,y)e xs+yt (36) 



where D(x, y) and E(x, y) are the Fourier transforms of D(s, i) and E(s, t) when n — > 0, 

D(x,y) = [ ^.±D(is,it)e- i ^+ t y\ E(x,y) = [ ±^E(is,it)e-^+^ , (37) 
J Zn Zn J Zn Zn 

and we expand all quantities at order n (i.e. A — Aq + nA\ + 0(n 2 ), E(x, y) = Eo(x,y) + nEi(x,y) + 0(n 2 )) as 
detailed in Appendix B. This allows one to solve the equations at the first two orders in n and, after some lengthy 
algebra, to compute Ai(g) (see Eq. B25) from which we finally get 



^q{9,H) = -[1 + b+ — }[S 3 {g)~g — ] (38) 

where 83(g) = J dxdyDo(x, y) ln[2 coshx]. A closed expression of m(g,H) (Eq. B27) may be also obtained from 
Eq. (22), but it contains more complicated integrals of Do(x,y) and it is in fact more convenient to compute the 
magnetization by numerically differentiating the expression (B25) of Ai(g) with respect to g. 

The only remaining task consists in solving numerically the coupled integral equations satisfied by Do(x,y) and 
Eg(x,y). This is outlined in Appendix C. 



IV. RESULTS FOR THE ONE-DIMENSIONAL CHAIN (z = 2) 

When z = 2 the Bethe lattice is just a collection of disconnected closed loops. In the thermodynamic limit, the 
dominant contributions to self-averaging quantities come from the loops of macroscopic size, and one thus expects 
the results to be the same as for the one-dimensional RFIM. Although no phase transitions occur in one dimension, 
frustration makes the problem nontrivial, and the relation between the shape of the hysteresis loop and the distribution 
of the metastable states can be fully elucidated. Note that the typical number of metastable states of the random- 
field Ising chain and their distribution in magnetization and energy were computed a few years ago by a numerical 
transfer-matrix method (l8l|. but the study was restricted to the case of zero external field: therefore, the relation with 
the T = hysteretic behavior was not investigated. Moreover, some results of Ref.^^| appear to violate the convexity 
inequality In < In AT, as will be shown in the following. 



A. Annealed complexity 

Let us first consider the total annealed complexity S™ aa: (iJ) as obtained from Eq. (26) with 5 = 0. Its dependence 
on H is illustrated in Fig. 1 for different values of the disorder strength R (hereafter we take J = 1 in the numerical 
applications) . It can be seen that the presence of the magnetic field reduces the total number of metastable states, in 
agreement with the expectation that a sufficiently large field imposes the existence of a unique state. As R decreases, 
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FIG. 1: Total annealed complexity T,J ax (H) for z = 2 and R = 1 (solid line), 0.5 (dotted line) and 0.25 (dashed line). 



the metastable states concentrate in the vicinity of H = and in the limit R — *■ they are completely destroyed by 
an infinitesimal field. From the solution of the equation t(v) = 1, the zero-field complexity is given by 



1 1 



<(H = 0)=]n[- + -\ 1 



(39) 



It may be noted that this quantity approaches In 1+ 2 V ^ 

1 + V5 



In- 



in agreement with a recent transfer matrix calculation 21 

0.188 when R — > 0, a value much smaller than the complexity of the pure ferromagnet, Y7^q{H = 0) - ... , 
0.481 [Il[l7|. This singular behavior of the zero-field complexity computed from Eqs. (25-28) in the limit R — > 
comes from the fact that marginal states, absent when R ^ 0, are not properly counted when R = 0: indeed, 
/(z) = [1 + erf(ar/iJv/2)]/2 -> 9(x) for x 7^ whereas /(0) = 1/2; to also include in our calculation the marginal 
states that exist in the pure system, one should instead set 8(0) = 1 as in Ref.0]. Actually, all metastable states are 
marginal in the Ising chain at zero field (they simply correspond to different positions of the walls separating domains 
of up and down spins) and one can directly calculate a "biased" complexity by associating to each domain wall an 

arbitrary weight p with < p < 1; the result is In 1+ ^ +4 ^ -, which for p = 1/2 identifies with the complexity found 
in the R — > limit of the present approach. 

In Fig. 2, we plot the projection of the annealed complexity with fixed magnetization E^(m, H) onto the H — m 
plane for R = 2. Only the region where the complexity is positive is shown: this is the region where the number of 
metastable states is exponentially large and thus comparable to the total number of states in the system; in the rest of 
the plane, E^(m, H) < 0, and the number of metastable states is equal to zero in the thermodynamic limit. Several 
curves in Fig. 2 are of special interest: i) the contour E^m, H) = 0, ii) the curve m^(i?) which is the locus of the 
maxima of Ea(w, H) and thus gives the magnetization of the metastable states dominating the (unrestricted) average 
number Af(H) in the thermodynamic limit (i.e. T,™ ax (H) = S^(myi(ff), H)), and iii) the T = saturation hysteresis 
loop calculated in Ref.0 from probabilistic arguments. We see on this example that the contour T*A(m, H) = 
significantly overestimates the size of the actual hysteresis loop (a similar observation was made in the case of the 
Sherrington-Kirkpatrick modeljlOj). The fact that the saturation loop lies well inside the contour E^m, H) = can 
be understood as follows. Each realization of the disorder has a different saturation loop and, as shown in Appendix 
A, there are no metastable states outside that loop. However, the observed saturation loop is obtained by averaging 
the magnetization over all disorder realizations, an averaging which is dominated by the typical realizations and as 
such is related to the quenched complexity Eq(?tt,, H) that from Jensen's inequality is always less than S^(m, H). So, 
metastable states do exist outside the (average) saturation loop, there may even exist an exponentially large number 
of them, but these states are non-typical and essentially non-observable. 
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FIG. 2: (Color online) Projection of £a(wi, H) onto the H — m plane for z — 2 and R = 2. The contours are separated by 
10 -2 (with Ea(0, 0) ~ 0.100). The red curve is the locus m,A{H) of the maxima of T,A(m, H). The blue curve is the T — 
saturation hysteresis loop calculated in Ref.jjj]. 




FIG. 3: Annealed (solid line) and quenched (dashed line) complexities in zero external field for z — 2 and R = 2 as a function 
of magnetization. 



B. Quenched complexity 

As discussed in section 3 and Appendix B, the quenched complexity and the magnetization are parametrized by g 
(see, e.g., Eq. (38)), with the maximum number of metastable states corresponding to g — 0. As g goes from to +oo 
(resp. -co), £q(<7, i?) decreases monotonically towards whereas m(g,H) increases (resp. decreases) monotonically, 
the surface Eq(to, H) having roughly the same overall shape as S^(m, H): this is illustrated in Fig. 3 for R = 2. 
Note however that the derivative <9E(m, H)/dm, i.e., the value of g, becomes infinite for nontrivial values of the 
magnetization for the quenched complexity whereas this only occurs for m = +1 or — 1 for the annealed complexity. 

The projection of T,Q(m,H) onto the H — m plane is shown in Fig. 4 for R = 2, together with m(g = Q,H), the 
typical magnetization of the metastable states [24|. The region where Eg(m, H) is positive is clearly thinner than 
the corresponding region for S^(m, H) shown in Fig. 2. In fact, as proven analytically in Appendix C, the contour 
T,Q(m, H) = obtained in the limit g —* ±oo corresponds exactly to the saturation loop calculated in Ref.Q. This 
means that the T — hysteresis loop is indeed the enveloppe of all the typical metastable states in the thermodynamic 
limit (i.e. the entropy density is strictly positive everywhere inside the loop). 

We show in Fig. 5 the total quenched complexity in zero field T,Q ax (H = 0), which is equal to Sq(to = 0, H = 0), 
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FIG. 4: (Color online) Projection of T,Q(m,H) onto the H — m plane for z = 2 and R — 2. The contours are separated by 
5.10 -3 {with Eq(0, 0) ~ 0.082). The contour ^Q(m,H) = (in blue) coincides with the saturation hysteresis loop calculated 
in Ref.jj], as shown in Appendix C. The dashed curve is the annealed contour Ea (m, H) = 0. The red curve is the typical 
magnetization of the metastable states (the locus of the maxima of the quenched complexity) . 




FIG. 5: Total quenched complexity in zero field T,Q ax (H = 0) for z = 2 as a function of disorder strength. 



as a function of R. It can be seen that the curve is very fiat near R = and that lim^^o Y,Q ax (H = 0) w 0.179 (see 

Appendix C), which differs from the limit of the annealed complexity liniR_>o T,™ ax (H = 0) = In 1+ g w 0.188 (the 
difference is well outside the error bars in the numerics), itself smaller than the total complexity of the pure system. As 
discussed above in commenting the latter result, those differences come from the singular features of the limit R — > 0: 
when R ^ 0, no metastable states are marginal, whereas when R = all metastable states are marginal. However, 
contrary to what happens for the annealed complexity in the limit R — *■ 0, the present procedure for calculating the 
quenched complexity does not reduce to a simple (and biased) way of counting the marginal states in the pure system. 
(Note also that Y^ ax {H = 0) sa 0.176 for R = 0.75, which is at odds with the value ^ ax {H = 0) » 0.25 that can be 
extracted from the Figure 5 of Ref.^^; this latter result is actually larger than the corresponding annealed complexity 
obtained from Eq. (39), YT2 ax {H = 0) « 0.186, which clearly violates Jensen's inequality and suggests some flaw in 
the calculations of Ref.[lq.) 

V. ANNEALED COMPLEXITY FOR z > 2 

For z > 2 one expects new phenomena to take place. In particular, the complexity may not be anymore a convex 
function of H and m and phase transitions may occur. 
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FIG. 6: Annealed complexity Ea(?ti, H) for z = 3 as a function of m for different values of H. Left panel: R — 2; from top to 
bottom = 0, 1.30, 2. Right panel: R = 1.2; from top to bottom H = 1.30, 1.34, 1.355, 1.372, 1.40 (the values of T, A (m = 1, H) 
not visible here are slightly negative). 

0.2 | , 1 , 1 , 1 , 1 




FIG. 7: Total annealed complexity Y.™ ax {H) for z = 3 and R = 2 (solid line), 1.5 (dotted line), and 1.2 (dashed line). Note 
the discontinuity in the slope of Y^ ax {H) at H « 1.355 for R = 1.2. 



Let us first consider the case z = 3. In Fig. 6 we show X^m, i?) as a function of the magnetization for two different 
values of R and several values of the external field (we only consider positive values of H because of the symmetry 
S(m, — H) = H(—m,H)). As can be seen in the left panel of the figure, the complexity is always a convex function 
of the magnetization when R is large enough. As a consequence, the total complexity E™ aa; (£f) = max m X^m, H) 
is a continuous function of H, as shown in Fig. 7. This is also the case for the corresponding magnetization mA(H). 
There are only two values of m for which XU(m, H) = (note that the two limiting values lim m ^±i XU(ra, H) = 
ln[(l + erf[(zJ ± H)/(R\/2)]) /2] are negative), and the overall surface XU(m, H) has roughly the same shape as for 
z = 2 (the projection onto the H — m plane resembles the one shown in Fig. 2). 

The behavior is more complicated at low disorder as illustrated by the curves in the right panel of Fig. 6. For R 
smaller than some critical value (R < 1.28) and in a certain range of H , XU(m, H) may indeed display two positive 
maxima at two distinct values of m. As a consequence, there is a discontinuity in dY72 ax {H)/dH at some value of 
the field (H w 1.355 for R = 1.2) as shown in Fig. 7 and, accordingly, there is a jump in rriA(H). Moreover, for 
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R even smaller, there is a small range of H for which Sa(w, H) — for four distinct values of the magnetization. 
This yields an overall surface E^(m, H) that has a more complicated shape than in the strong-disorder regime, as 
illustrated by the projection onto the H — m plane shown in Fig. 8. In particular, one can observe in the right 
pannel the "S" shape of the contour E^m, H) = in the range 1.3 < H < 1.4. Note that the actual hysteresis loop 
calculated in Ref.JJ is again inside the contour E^m, H) = as it must be. The fact that drriA(H) / dH is negative 
in a certain range of H is surprising and can only be explained by the existence of metastable states that have a quite 
non-typical magnetization. One must recall, however, that rriA{H) is not an observable quantity. For the same reason, 
the non-monotonic behavior associated with the contour E^(m, H) = bears no connection with the behavior of the 
actual hysteresis loop: this latter does not display any finite jump nor phase transition for z = 30. 




FIG. 8: (Color on line) Projection of Y*A(m,H) onto the H — m plane for z = 3 and R — 1.2. In the left panel, the contours 
are separated by 2.10 -2 (with Ea(0, 0) ~ 0.171). The red curve is the locus iiia(H) of the maxima of Ea(ti, H). The blue 
curve is the T = hysteresis loop calculated in Ref. 4]. The right pannel is a close-up of the figure in the range 1.3 < H < 1.4, 
with a spacing 5.10~ 4 between the contours. 

For z = 4, the situation changes for the hysteresis: a phase transition and a jump discontinuity appear for R < 
1.780. E>t(m, if), on the other hand, has the same type of behavior as for z = 3. In particular, for R < 2.25 and 
in a certain range of the external field if, it has two positive maxima as a function of m and four zeros, like in the 
left panel of Fig. 6. The projection of E^m, if) onto the if — m plane is shown in Fig. 9 for R = 1.5. mA(H) has 
again a strange behavior, first decreasing (resp. increasing) as H is increased (resp. decreased) and then jumping to 
a value which is very close to 1 (resp. —1). This is probably the generic behavior for z > 4. 



VI. DISCUSSION 



One important conclusion of the computations presented above is that the annealed approximation does not properly 
describe the actual distribution of the metastable states of the RFIM in the field-magnetization diagram (i.e., the 
distribution of metastable states associated with typical realizations of the disorder) . The calculations are done here 
on a Bethe lattice, but this conclusion is probably true for any lattice. Predictions can even be qualitatively wrong, 
as illustrated by the fact that a phase transition corresponding to a change of convexity of the annealed complexity 
occurs for z = 3 whereas the typical behavior changes for z > 4 only. On the other hand, the quenched calculation 
for z — 2 shows that the entropy density associated to the metastable states is strictly positive everywhere inside the 
saturation hysteresis loop. It is very likely to be always the case in the strong-disorder regime where the saturation 
loop is continuous. The situation is less clear for z > 4 in the weak-disorder regime where the saturation loop displays 
a jump discontinuity in the magnetization. Since we have not so far been able to carry out the quenched calculation 
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H 



FIG. 9: (Color on line) Projection of TiA{m,H) onto the H — m plane for z = 4 and R = 1.5. The contours are separated by 
1.5.10 -2 (with £a(0, 0) ftf 0.159). The red curve is the locus mA(H) of the maxima of £a(?ti, H). The blue curve is the T = 
hysteresis loop calculated in Ref.Q. 



m 




H Hi H2 

FIG. 10: (Color on line) Speculative distribution of the metastable and H-states in the H — m plane for z > 4 in the weak- 
disorder regime where the hysteresis loop has a jump at H = ±H%. The number of metastable states is exponentially large 
in system size in the shaded region and zero outside. Dots indicate the region where the number of H-states is exponentially 
large. The red curve is the typical magnetization of the metastable states. 

to the end, we can only speculate about the actual behavior of the quenched complexity T,Q(m, H). A possible 
scenario is depicted in Fig. 10. The rationale for this scenario is very intuitive: we simply assume that the jump in 
the magnetization curve at some field H2(R) (or, symmetrically, at —H2(R)) occurs because, with probability 1 in 
the thermodynamic limit, there are no metastable states with a low magnetization for H > Hi- The "S" shape of 
the curve Xq(?ti, H) = is speculative but it is tempting to associate this contour with the real but "unphysical" 
roots of the self-consistent polynomial equation for the fixed-point probability P* derived in Ref.[4j to describe the 
saturation hysteresis loop. Let us recall that for z — 4 this equation has three real solutions in the range H\ < H < H2 
with two of them merging and becoming complex at H — Hi and H — Hi ■ Of course, one cannot discard a more 
complicated scenario to occur, in particular when there are more than three roots for z > 4. In principle, there is 
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also no objection for the existence of a jump discontinuity in the typical magnetization, i.e., the magnetization of the 
states corresponding to max m Eq(to, H) (which is not the case considered in Fig. 10). One could indeed imagine 
a mechanism similar to the one displayed in Fig. 6 with Eg(m, H) having two maxima as a function of m in a 
certain range of H: a "first-order" transition could then occur when the two maxima have the same height. These 
are interesting open questions which remain to be resolved before the problem of the distribution of metastable states 
of the RFIM on the Bethe lattice can be considered as fully understood. It would be of course even more interesting 
to know the actual behavior on a 3-d lattice, but this can only be done numerically. 

Finally, let us remark that very little is known about a special category of metastable states, the so-called "in- 
states" |26j , which are the spin configurations that can be reached by a field history starting from one of the saturated 
states (these configurations are characterized by the applied sequence of reversal fields {Hi}). The Zf-states are thus 
experimentally observable by recording various hysteresis loops and subloops. The most easily accessible ones, which 
we have discussed in this paper, are those associated with the saturation hysteresis loop. The ff-states are believed to 
be low-energy configurations, and it has been recently suggested that their number also increases exponentially with 
the system sizej23|- An analytical check of this result would be much welcome, even in one dimension. As was pointed 
out in introduction, in the weak-disorder regime, for z > 4, a region around the origin in the (H — m) plane becomes 
inaccessible to any field history, which forbids a demagnetization process to be completed |l2j. The corresponding 
distribution of the if -states is as illustrated in Fig. 10. In particular, the probability of finding iJ-states in the central 
part of the loop (the region of the (H — m) plane without dots) is zero in the thermodynamic limit. This is related 
to the fact that the states on the concave parts of the envelope of the shaded region are not reachable by any field 
history because of the jump in the major hysteresis 1oop|27|. 



We acknowledge useful discussions with V. Basso, S. Zapperi, and E. Vives. The Laboratoire de Physique Theorique 
des Liquides is the UMR 7600 of the CNRS. 



In this Appendix, we show that the probability to find metastable states outside the saturation hysteresis loop 
vanishes in the thermodynamic limit when the interactions are ferromagnetic. This results from the fact that for 
any given realization of the random fields, there can be no metastable states outside the corresponding saturation 
hysteresis loop. Indeed, the properties of the T = single-spin-flip dynamics are such that for any given value of the 
magnetic field H, all sequences starting from the fully polarized configuration (with all spins up or all spins down), 
sequences in which unstable spins are flipped, lead to the same metastable state on the upper (resp. lower) branch 
of the hysteresis looppj. This is true in particular for a parallel updating of all unstable spins. Let us denote C a (H), 
a = l,..,,p, the p configurations obtained after each updating steps along this process with Cq(H) being the initial 
fully polarized state with, say, all spins up and C P (H) — C+(H) the final state on the upper branch of the hysteresis 
loop. The magnetization of course decreases as a increases. Consider now a metastable configuration C(TL). Because 
of the ferromagnetic nature of the couplings, the unstable spins that flip from Cq to C\ must be down in C(H); the 
same is true for the unstable spins that flip between C\ and C2 and the reasoning can be repeated so that all down 
spins in C P (H) = C+{H) must also be down in C(TL). The magnetization of C(TL) is thus less than that of C+(H). 
This is a form of the "no-passing rule" [I| . The argument also applies when starting from the fully polarized state 
with all spins down, which completes the demonstration. As the property is true for all realizations of the disorder, 
this applies in particular to the "typical" realizations whose number dominates the whole distribution, so that one 
expects no "typical" metastable states outside the (average) saturation hysteresis loop. 
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APPENDIX A 



APPENDIX B 



In this Appendix, we detail the calculation of the quenched complexity for z 
Injecting (35) into Eq. (34), we get the following equations, 



2. 




1 - a 
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[A + Be~ 9n + C + £>(-!,-!) + ^ E(n, 
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(Bla) 
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— C = \^e~ gn [A + Be~ 3n + C + E(l, 1) + ^ D{n, 1)] + ^±£ e -" n C (Bib) 

-B=^C+ b -±^A (Blc) 



and 



D(tr,T) = ?—^-e as {Be gs + D{a,<r) + ^E(Ti,a) 

2 Ti 

x J] [<*K; -l)^i a ; -r a ) + *W;r»)i(/; a")] } (B2a) 



— £(<r, r) = £±£ e »«{Be»» + £(<t, <x) + ]T D(n, <r) 

x [] [S(a a ; l)6tf; -r a ) + 5(^;r a )5(T a ; a a )] } . (B2b) 

a 

After some algebra we also obtain from Eq. (21) the additional equation 

{A + Cf + 2B{Ae~ gn + Ce gn ) + B 2 (2 cosh 2g) n + 2BJ2 e 9 " [D(<r, tr) + E(tr, a)] 

(T 

+ 2A[D(-1, -1) + E ( T > - 1 )] + 2C[S(1, 1) + 2 ^> !)] 

T T 

+ ^[£> 2 (<r, <r) + E 2 (ct, <r)] + 2j2 D(a, t)E(t, ct) ]J[6(a a ; l)S(r a ; -1) + S(a a ; r a )} = 1 (B3) 



er.r 



Eqs. (B1-B3) constitute a set of six coupled equations that must be solved at the first two orders in n. As already 
mentioned after Eq. (35), it is important to note that D(cr, t) and E(cr, r) are functions of two variables only. Indeed, 
the term ]J a [S(a a ; -l))<5(r a ; l) + <5(r a ; cr a )] in Eq. (36) imposes that a a = 1 and r a = -1 cannot hold simultaneously. 
This in turns implies that u = n — \s — t\. 

We now expand both sides of Eqs. (B1-B3) at order n. Let us first consider objects like D(—l, —1), ^2 T D(n, 1), 
etc... in Eqs. (Bl). Using the integral representations (37) of D((t,t) and E(ct 7 t) 7 wc find 

D(-l, -1) = J dxdyD(x,y)e- n( - x+ ^ 

= D Q + n{D 1 - Si - S 2 ) + 0(n 2 ) (B4) 



5^D(ti,1)= f dxdyD(x,y)e n yJ2 

Ti J Ti 



e *EaTl 



= D Q + n(D 1 + S 2 + S 3 ) + 0(n 2 ) (B5) 

where Do — J dxdyDo(x,y), D\ = J dxdyDi(x,y), Si — J dxdyxDo(x,y), S2 = J dxdyyDo(x,y), and 53 = 
/ dx dy Do(x,y)\n[2 cosh x]. (Similar expressions are obtained for £7(1,1) and 5Z r E(n, -1), with Eo, Ei, S[, S' 2 
and S' 3 replacing respectively Do, D\, — Si, — S 2 and S3.) 
In the same way, one has in Eq. (B3) 

^e«'[D(a,a)+E{a,a)]= [ dxdy[D(x,y) + E(x,y)}J2^ i9+X+v) 

= Do + E + n(Di + Ei + S 4 ) + 0{n 2 ) (B6) 
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e s(x+y+x' +y') 



Y^[D 2 {cr 7 cr) + E 2 (cr } cr)} = f dxdydx' dy'[D(x, y)D(x', y') + E(x, y)E(x', y')] ^ . 
cr •* cr 

= D 2 + El + n(2D D 1 + 2E E 1 + S 5 ) + 0{n 2 ) (B7) 



J2 D{a, t)E(t 7 cr) l[[S(a a ; -l)5( T a ; 1) + S(a a ; r a )} 

T a 

j dxdydx 1 dy'D(x, y)E(x', y') ^ e<"+»')+t(y+x') Y[[S(a a ; -l)S(r a ; 1) + S(a a ; r a )] 

■' i-r t- 



CT,T 



ct,t 

2 



= D E + n(A)£i + DxE + S 6 ) + 0(n 2 ) (B8) 

where S4, = J dxdy[Do(x,y) + Eq(x, y)} ln[2 cosh(ir + y + g)], S5 = J dxdydx 1 dy'[Do(x,y)Do(x' ,y') + 
E (x,y)E Q (x', y')]ln[2cosh(x + y + x' + y')] and S 6 = J dxdydx' dy' D (x, y)E {x' ,y')\n[ e - x+v+x '- v ' + 2cosh(x + 
y + x' + y')]. 

The treatment of Eqs. (B2) is slightly more involved. Consider first the term T = 

e Ba ^ Tl E(T 1 ,a)lla[ 5 ( aa \- 1 ) S ( T ii- Ta ) + 5(T?;T a )5(T a ;a a )] in Eq. (B2a). Performing the sum over t x 
gives T = J dxdyE(x, y) \[ a f a * T * (x, y) with 

f a « T a(x,y) =e-"°-»-'%Vl) + e (l+s+9) '°%V a ) • (B9) 

Following Ref.^Jj we now use the identity 

JJ/ ff . T .=expEln/ ff . T .] =exp[^^ 6(a; a Q )<5(r; r a ) ln/ CTT ] (BIO) 



a cr,r 



with 

7 (s + u + < + n) for cr = t = 1 



^ v 1 ±(— s — u + i + n) for cr = — l,r=l 



i(— s + u — t + n) for a = r = — 1 . 
Remembering that the combination cr a = 1, r a = — 1 is forbidden and that u = n + s — t, we then obtain 

JJ/ CT a T a(x,y) = exp[^t^ln/ ++ (a;,y) + ^— - ln/_ + (a:,y) - t — — ln/__(i, y)] (Bll) 

a 

with f++{x,y) = e x+ y+9, f_+(x,y) = e -^ x+v+ ^ and = e * + e-^+f+s). Hence 

T = / dxdyE(x, y) exp[s(a; + y + .g) - - ln(l + e 21 ') + - ln(l + e 2 *)] . (B12) 
Using the same procedure, we get 

e 9S D(er,er) = J dxdyD(x,y) exp[s(x + y + g)} . (B13) 

Finally Fourier transforming both sides of Eq. (B2a) with respect to s and t, we obtain 

e A - , , a - 6 



-£>(x, v) = — {BS(x - 2g)S(y) + S(y) J dx> ' D{x' , x-x' - g) 

+ J dx'dy'E{x',y')5(x- x' -y' - g)5(y+^\n[l + e 2x '})e% ln(1+e2x '^} (B14) 
In the same way the Fourier transform of Eq. (B2b) gives 

?-E(x,y) = ^{B5(x-2g)8(y) + 5(y) J dx'E(x', x - x' - g) 

+ f dx'dy'D(x',y')5(x-x' -y' -g)S{y-lln[l + e- 2x '])e^ Mi-e 2 *')} ( B15 ) 
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Collecting all terms at order n°, we obtain the six coupled equations 



e 



Ao 



1 — CL r . „ „ „ „ , a — b 



and 



= -g-IA) + S + Co + A) + Eo] + —A (B16a) 



— Co = i-^o + Bo + 0) + Do + £b] + ^C (B16b) 

<y 2 2 



—So = ^Co + (B16c) 
7 2 2 

= ^{.Bo<5(:r-2 3 )%)+%) / dx' D (x' , x - x' - g) 

+ J dx'dy'E Q {x', y')8{x - x' - y' - g)S(y + X - ln[l + e 2x '})} (B17a) 

E (x,y) = ^{B S(x-2g)6(y)+S(y) J dx'E (x', x - x' - g) 

+ Jdx'dy'D (x',y')6(x-x'-y'-g)6(y-±ln[l + e- 2x '})} (B17b) 

(A + B + C + D0 + E0) 2 = 1 . (B18) 



Eqs. (B16) and (B18) may be decoupled from Eqs. (B17) by performing the integration of (B17) over x and y which 
readily gives 

a — b b + c ,„ 

2 — (a + c) 2 — (a + c) 

The only physical solution of Eqs. (B16), (B18) and (B19) is then 

A = ln 7 (B20a) 



A °=2^b) ^ 

_ (a-b)(l-c) (b + c)(l-a) 
B ° 4-2(6 + c) + 4-2(a-6) (B2 ° C) 

Co = . (B20d) 

Note that Eq. (B20) is necessary for the n — > limit to exist when z = 2 (as pointed out after Eq. (30)). This solution 
corresponds to the choice A Q + B Q + C a + D Q + E Q = 1 in Eq. (B18) which implies that lim„^ J2er r c ( cr ' T ) = \/2/7j 
a result that could have been also obtained from Eq. (21) by naively (but incorrectly) using a Fourier integral 
representation of c(er, r). 

It remains to solve the coupled integral equations (B17). This can be done by writting formally Do(x,y) and 
E (x, y) as two infinite series of delta functions. The locations of the peaks are then obtained recursively as discussed 
in Appendix C. 
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We now achieve the calculation of the quenched complexity by solving Eqs. (B1-B3) at first order in n. Using Eqs. 
(B4-B8), (B14-B15), (B17-B18), and (B20a), we obtain, after some manipulations, 

[Ai + l ~^g]A + = 1^ [ - g(C + 2B + D Q + E ) + C 1 +B 1 +D 1 +E1 

-Sl-S2 + S' 3 - S' 2 ] (B21a) 

[Ar - ^g]C + i^Ci = ^ [g(A + 2B Q + D + E ) +A 1 +B 1 +D 1 +E 1 

+ S[+S' 2 + S 3 + S 2 ] (B21b) 



ArBo + B 1 = + h ~^M (B21c) 

D^x, y) = -(y + Ai)D (x, y) + ^{B^x - 2g)5{y) + S(y) J dx'D^x', x-x' - g) 

+ j ' dx'dy'Eiix^y'Wx-x' - y' - g)5{y +^\n[l + e 2x '])} (B22a) 

E 1 {x,y) = {y-K 1 )E {x,y)+ h -^{B 1 8{x-2g)5{y) + 5{y) J dx'E x (x',x - x' - g) 

+ J dx'dy'D^x', y')5(x - x' - y' - g)S(y - \ ln[l + e - 2x '})} (B22b) 

A 1 +B 1 + C 1 +D 1 +E 1 + g(C - A )B + ^ ln[2cosh(2 3 )] + (C - A )(S 2 + S' 2 ) 
+ C (S[ + S 3 ) + A (S' 3 - + B S 4 + ^ + S 6 = . (B23) 
Integrating Eqs. (B22) over x and y we obtain 

D 1 = a _~ h _ (Bi - AxC + S' 2 ) - 2 ~ b ~ C (S 2 + AM (B24a) 

A CI C A CI C 

Et = b + C {B 1 - A.Do - S 2 ) - l~ a + b (S 2 - A.Do) (B24b) 

2 — a — c 2 — a — c 

which again allows to decouple Eqs. (B21) and (B23) from (B22). The solution of the system of linear equations 
(B21), (B24), and (B23) then yields 

Ai(.9) = ^{g[(b-c)A + (a + b)C + (a-c)(2B + D + E Q )] - S^l - a) + S[(l - c) 

+ S 2 (a -c-2) + S' 2 (a-c + 2) + S 3 (l - c) + S' 3 (l - a)} (B25) 

which, remarkably, does not depend on the integrals S4, S5 and 5*6. In fact, it turns out that Ai(g) can be written in 
terms of a single integral of D {x, y) or E (x, y) (for instance S3). Indeed, the analysis of Eqs. (B22) shows that 

Si = — ■ -g , S 1 = -Si B26a 

(1 -p- q)(l -p- q + pq) p 

S 2 = + S 3 ) , S' 2 = |(S 3 - Si) (B26b) 
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S' 3 = ^S 3 . (B26c) 

with p = (a-b)/2 and q = (b + c)/2. Inserting (B26) in Eq. (B25), using Eqs. (B19) and (B20) and the relation 
EqG?, H) = Ai(g) - gdA^/dg, we finally obtain Eq. (38). 

Note that the magnetisation m(g) in the limit n — > only depends on the solution of the saddle-point equations at 
lowest order. Inserting Eq. (35) in Eq. (22) gives after some algebra 

m(g) = Cl - A% + 2B (C - A ) + 2{C Q E - A Q D Q ) + B% tanh 2g 

+ A n J dxdyE a (x, y) [tanh a; — 1] + C J dxdyD a (x, y)[tanhx + 1] 

+ 2B J dxdy[D (x,y) + E Q (x, y)} tanh(x + y + g) 

+ J dxdydx'dy'[D (x, y)D (x' , y') + E (x, y)E {x' , y 1 )} tanh(.x + y + x' + y 1 ) 

/ x+y+x'+y' _ -(x+y+x'+y 1 ) 



APPENDIX C 

The solution of the coupled integral equations (B17) cannot be obtained in the form of closed expressions for Dq(x, y) 
and E (x,y). However, this type of self-consistent equations is common in one-dimensional disordered systems and 
it is easy to see that Do(x,y) and Eo(x,y) may be represented as two infinite series of delta functions which are 
conveniently written as 

oo i 

D (x,y)/B = pj2Y,p J( i t ~ l£) v( x >y) ( Cla ) 

CO i 

E (x, y)/B = qJ2J2 tfP^Eijix, y) (Clb) 

with p = (a - b)/2, q = (b + c)/2 and 

0) 

Dij(x,y) = J2 S ( X - 4 3) (9))S(y y^Hg)) (C2a) 
k=i 

G) 

E ij (x,y) = J2s{x-x^\g))5{ y - y , ^\g)) . (C2b) 
k=i 

L>ij(x,y) and Eij(x,y) are then computed iteratively, the coordinates (x^\y^) and {x^\ y '^) of the delta peaks 
at step i being obtained from the coordinates at step with Doq(x, y) = Eqo(x, y) = 6(x — 2g)S(y). For instance, 

the first terms in the expansion of Do(x,y) are 



Dio(x,y) 


= S(x 


-3g)6(y +^ln(l + e^)) 


Du(x,y) 


= S(x 


- 3g)S(y) 


D 2 o(x,y) 


= S(x 


-Ag)8{y+ 1 -\n{l + e^)) 


D 2 i(x,y) 


= S(x 


-Ag+ l -\n{l + e^))5(y) 


D 2 2(x,y) 


= S(x 


- 4<?)<%) 


etc... 







l -Hl + e-^))5(y+ l -Hl + e^)) 



(C3) 
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Note that p and q depend on H and R whereas the dependence on g is only contained in the position of the delta 
functions. With Eq. (Bla), the integral £3 also becomes an infinite series, 

00 i G) 

S a (g) -SopEE E M 2 cosh x k 3) (»)] ( C4 ) 

1=0 j=0 k=l 

The number of delta functions grows exponentially during the iteration procedure and the precision on S3 depends 
on the number / of terms that are included in the infinite sums. However, the cumulative weight of the delta peaks is 
finite (D /B = p/(l — p — q) and Eq/Bq = q/(l — p — q)), and one can take as many terms as needed so to satisfy a 
given accuracy defined by e = 1 — (1 — p — q) J2i=i £}=o (j)^?*^ ( see > e -g-> Ref.pEj for a similar calculation). When 
H = 0, however, one has p = q = (1/2) erf(v / 2J/-R) and the convergence of the series becomes progressively poorer 
as the ratio J/R increases and p, q — > 1/2. In this case, it is convenient to write Sa(g) as 

00 

s 3 (g) = Bojz^-ssig) + A>pZy [4%) - '3(3)] (C5) 

where s^\g) — J2]=o Sfe=i m P cosh (g)] and s%°(g) — limj—nx, s$\g). For g = 0, for instance, one finds numeri- 
cally sf(0) re 0.716 and since B ~ 1 - 2p when R->0, this yields from Eq. (38) lim fl ^ Zq cix (H = 0) = s§°(0)/4 ps 
0.179. 

When g — > ±00, the expressions of Do(x,y) and Eo(x,y) simplify considerably as the positions of the delta peaks 
concentrate on a few lines in the x — y plane. Specifically for g — > +00, one obtains 

A>(x, y)/fl = V[p"- l (1 ~ P)(1 ~ 9) ^ - n 5 )%) + - l)5(x - ng)6(y - (1 - n) 5 )] (C6a) 
^2 1-P"« 



with 



Replacing in Eq. (B27) gives, after some algebra, 



00 

Bb(s, y)/B = E E(n)S(x - ng)S(y) (C6b) 

n=2 



* ( „,_ l £frtaz*i^ +1 ^rL^ 1 ( C7) 



m(.g -> +00) = 1 - 2A (1 + B + D ) 

= 1 - 2p(-H) (C8) 



with 

, (a- l)(a + ac-4 + 6 2 -2&-c) 
P[ - H) = (2-a + b) 2 

1 + PQP2 ~ P\ fnQ , 

= P°TT I 12 C9 ) 

[1 + po -pi] 1 

where p n (—H) = f^ri-nt+H P(h)dh is the probability (introduced in Ref.Q) that a spin with n neighbors up (n — 
0, 1,2) is up at the field -H: p = (l-a)/2,pi = (l-6)/2,p 2 = (l + c)/2. The magnetization given by Eqs. (C8-C9) 
is exactly that calculated in Ref.0 for the upper branch of the hysteresis loop. By symmetry, the magnetization on 
the lower branch of the loop is given by m = 2p(H) — 1 and corresponds to the limit g — > —00. 

It remains to calculate Ai(g,H) and T,Q(g,H) = Ai(g,H) — gm = h\{g,H) — gdAi(g, H)/dg when g — > ±00. 
Inserting Eqs. (C6-C7) in (B25) gives, after some algebra, 

. . TT . (a 2 c- a - 2ac + ab 2 - 2b 2 - 26 + c) 

Ax(9 - +00, H) ~ (2 _ a + 6)2 " > (C10) 
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from which one readily finds that £q(<7 — > +00, H) — 0. By symmetry, Eq also vanishes along the lower branch of 
the hysteresis loop. To illustrate how the cancellation of terms occurs in £q(<7, H) in the limit g — > ±00, let us give 
the first terms in the large- R expansion of Ai(g, H) and £q(<7, H): 

Al(g,H) = 9\[\^ + ^i 2 9H + Jln(2 cosh2«?)] - _L_y -[^ff 3 + 6ngHJ 2 - 24gHJ 2 

J 4 

+ 24J 3 ln(2 cosh 2g) - 24J 3 ln(2 cosh3g)] + 0{ — ) (Cll) 

R 

2J 2 4J 3 [2 

T, Q {g,H) = ^[ln(2cosh25) - 2#tanh25] + ? A/-[2.gtanh2g - ln(2cosh2.g) 

ttR 2 ' nR A V 7r 

j4 

+ ln(2cosh3g) - 35tanh3p] + 0( — ) (C12) 

R 

Note that Eq(<7, i?) does not depend on H at the first two orders of the expansion. 
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